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Abstract 

Background: Comparing metabolic profiles under different biological perturbations has become a powerful 
approach to investigating the functioning of cells. The profiles can be taken as single snapshots of a system, but 
more information is gained if they are measured longitudinally over time. The results are short time series 
consisting of relatively sparse data that cannot be analyzed effectively with standard time series techniques, such 
as autocorrelation and frequency domain methods. In this work, we study longitudinal time series profiles of 
glucose consumption in the yeast Saccharomyces cerevisiae under different temperatures and preconditioning 
regimens, which we obtained with methods of in vivo nuclear magnetic resonance (NMR) spectroscopy. For the 
statistical analysis we first fit several nonlinear mixed effect regression models to the longitudinal profiles and then 
used an ANOVA likelihood ratio method in order to test for significant differences between the profiles. 

Results: The proposed methods are capable of distinguishing metabolic time trends resulting from different 
treatments and associate significance levels to these differences. Among several nonlinear mixed-effects regression 
models tested, a three-parameter logistic function represents the data with highest accuracy. ANOVA and 
likelihood ratio tests suggest that there are significant differences between the glucose consumption rate profiles 
for cells that had been-or had not been-preconditioned by heat during growth. Furthermore, pair-wise f-tests 
reveal significant differences in the longitudinal profiles for glucose consumption rates between optimal conditions 
and heat stress, optimal and recovery conditions, and heat stress and recovery conditions (p-values <0.0001). 

Conclusion: We have developed a nonlinear mixed effects model that is appropriate for the analysis of sparse 
metabolic and physiological time profiles. The model permits sound statistical inference procedures, based on 
ANOVA likelihood ratio tests, for testing the significance of differences between short time course data under 
different biological perturbations. 




Background 

Innovations in molecular biology, miniaturization, and 
robotics have led to genetic, proteomic and metabolo- 
mic data in quantities never seen before. Microarrays 
exhibit the expression state of thousands of genes in a 
single experiment, 2-D gels and other proteomic meth- 
ods identify subtle changes in the protein profile of an 
organism under altered conditions, while techniques of 
mass spectrometry are capable of quantifying hundreds 
of metabolites in one spectrogram. Most of these new 
tools have been utilized by biologists to obtain snapshots 
of the state of a cell population, and these snapshots 
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have been compared to controls in order to gain 
insights into the responses of cells to various stimuli or 
perturbations. Thus, a typical experiment might com- 
pare the gene expression in a cancer cell to the corre- 
sponding expression profile in normal, healthy cells. 

Comparisons between molecular profiles of perturbed 
cells and controls have vastly extended our understand- 
ing of the healthy and pathological functioning of cells. 
Their drawback is the singular nature of measurements 
at only one time point, which limits insights into the 
scope of a response and leaves doubt whether the best 
time point had been selected for data collection. 
Responding to this disadvantage of snapshot compari- 
sons, recent technical advances in molecular biology 
have rendered it possible, as well as practically and eco- 
nomically feasible, to execute high-throughput 
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measurements in time series. For instance, many data 
are now available characterizing the genomic responses 
of microorganisms to a variety of stimuli over a period 
of time [1]. Similarly, non-invasive nuclear magnetic 
resonance (NMR) measurements can characterize the 
time development of metabolites in intact cells over a 
span of hours, with data taken every minute or in some 
cases even every 10 to 20 seconds [2]. These time series 
data contain enormous amounts of information, because 
they are systemic as well as dynamic and capture the 
cellular response in a more comprehensive fashion than 
measurements at single time points can achieve. 
Furthermore, time series data constrain each other in 
some sense, for instance by facilitating the identification 
and interpretation of outliers and measurement errors. 

Biological time series data contain more information 
than snapshot data but they cannot necessarily be evalu- 
ated with the typical methods of time series analysis, 
such as autocorrelation and frequency domain techni- 
ques, because these were often developed for much 
longer and denser time series, as they are found, for 
instance, in weather recordings. In comparison to such 
data, presently available biological data are very sparse. 
Indeed, many time series consist of maybe a dozen 
quantities, measured at a handful of time points, and 
replicates are usually scarce. This paucity of data, 
accompanied by minimal redundancy, immediately cre- 
ates statistical challenges that have not been addressed 
in a systematic fashion. For instance, in the analysis of 
microarray time series [3], genes are often merely clus- 
tered visually into time trends, such as monotone 
increase/decrease, S-shaped increase/decrease, or initial 
increase/decrease followed by a gradual return to the 
initial state. For metabolic time series, even fewer meth- 
ods are available, and it is not often possible to state 
with a sufficient level of confidence and objectivity 
whether two time series are significantly different. 

We encountered this issue in the analysis of in vivo 
NMR time series measurements in the baker's yeast Sac- 
charomyces cerevisiae [4] and the lactic acid bacterium 
Lactococcus lactis [5], Investigating these systems under 
many different conditions, we found it difficult to decide 
how much of the differences between two time courses 
were due to normal variation between cell cultures and 
experimental settings and how much was pointing 
toward true differences in the cells' response dynamics. 
We therefore decided to develop sound statistical meth- 
ods for assessing the significance of differences between 
biological time series. 

Specifically, we analyzed the following situation. We 
grew yeast cells under optimal conditions to a particular 
population density and measured the uptake of their 
favorite substrate, glucose, over time. We compared 
these time-dependent uptake profiles with those of the 



same cells exposed to persistent heat, which in yeast 
evokes a systemic heat stress response. We then 
returned the cells to optimal temperature and measured 
the uptake dynamics under recovery conditions. The 
simple question we asked was whether the uptake char- 
acteristics are significantly different among the three 
situations of ambient optimal temperature, heat, and 
recovery. In addition, we were interested in answering 
the question of whether preconditioning the cells with 
heat during growth would affect their responses to heat 
later in life. 

While the biological implications are of relevance in 
themselves, the more important focus here is on more 
or less generic analytical tools for scarce time series in 
biology. Thus, we propose in the following the adapta- 
tion of statistical techniques to the significance analysis 
of differences between relatively sparse biological time 
series data. 

Methods 

Biological Methods 

Strain and growth conditions 

Saccharomyces cerevisiae strain JK93da was kindly sup- 
plied by Dr. Ashley Cowart, Medical University of South 
Carolina. Cells were grown up to an OD 60 o (cell density) 
of approximately 2 in GO minimal medium [6] supple- 
mented with 0.2% yeast extract, in a 5-liter fermentor 
(100 RPM, pH 6.5 under continuous flushing with air, 
pC>2 > 80%). We compared two sets of preconditions. 
Under control conditions, the cells were grown at an 
optimal temperature of 30°C throughout, whereas cul- 
tures preconditioned by heat stress were prepared as fol- 
lows: cells were initially grown under control conditions 
(30°C) up to an OD 60 o of approximately 1.3; at this 
point the temperature was increased to 39°C (heat 
stress) and the cells were allowed to grow for further 40 
min to a final OD 60 o of about 2. The actual experiments 
followed afterwards. 

Online measurements of glucose uptake by carbon-13 NMR 
( 13 C-NMR) 

Cells were harvested, centrifuged (10 min, 8,600 x g, 4° 
C), washed (5 mM potassium phosphate buffer (KPi), 
pH 6.5), and re-suspended in 50 mM KPj buffer con- 
taining 6% (v/v) H 2 0 and antifoam agent. The cell sus- 
pension (typically 60 mg dry weight/mL) was transferred 
to an 80-mL fermentor, maintained at 30°C, and con- 
nected to a 10-mm NMR tube by a circulating system 
[7]. The culture was pumped through the NMR tube at 
a rate of 30 mL/min and the pH was maintained at 6.5 
by automatic addition of NaOH or HC1. At time zero, 
[1- 13 C] glucose (65 mM final concentration) was added 
and the glucose concentration monitored by 13 C-NMR 
until substrate depletion, with a time resolution of 30 
seconds. The glucose was typically consumed within 10 
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minutes, and the cells began to starve for about 20 min- 
utes. Afterwards, the temperature was raised to 39°C, 
and a second pulse of glucose was added. Once glucose 
was exhausted again, the temperature was set back to 
30°C and a third pulse of glucose supplied. The 1 C- 
NMR time series were acquired on a Bruker DRX500 
spectrometer using a quadruple-nucleus probe head as 
previously described [6]. 
Experimental Data 

Six glucose uptake experiments were performed. They 
were essentially the same, but differed slightly in the 
amount of biomass at harvesting and, more importantly, 
in the temperature of the medium, which was initially 
optimal (30°C), then elevated to heat stress (39°C), and 
again reverted back to optimal (30°C) during recovery 
after heat stress. In the following, the experiments will 
be denoted as NG1 (normal growth; 83.85 mg/ml), NG2 
(normal growth; 57.05 mg/ml), NG3 (normal growth; 
60.3 mg/ml), NG4 (normal growth; 61.55 mg/ml), PCI 
(preconditioned; 73.45 mg/ml), and PC2 (precondi- 
tioned; 42.8 mg/ml), where the values in parentheses 
correspond to the biomass concentration used in each 
experiment. 

The cultures for experiments NG2, NG3, and NG4 
were grown to the target optical density (OD) of 2.0 and 
hence the respective biomasses are similar. Cells for 
NG1 were grown under the same conditions but har- 
vested later (OD = 2.8), thereby starting the actual 
experiments with a larger population size. As a conse- 
quence, the substrate was used up faster in this case. In 
order to check whether there are differences in future 
metabolic profiles of the cells which are grown under 
stressful conditions some of the cells (experiments PCI 
and PC2) were heat stressed during the final 40 minutes 
of growth prior to the in vivo NMR experiments (see 
also [4]). In all experiments, the glucose concentrations 
were recorded at equal time intervals of 30 seconds, 
starting roughly at times 3, 33 and 63 minutes for opti- 
mal, heat stress and recovery conditions, respectively. 
Figure 1 shows the original glucose consumption curves, 
stating when glucose was added (arrows), and Figure 2 
the corresponding rates of change in glucose concentra- 
tion for each experiment for different time points. As 
explained below, the time profiles were centered for the 
statistical analysis. 

Statistical Methods 
General objective 

The main objective of this study is to develop statistical 
tools that allow a biologist to decide whether differences 
in similar time courses (here, the change of glucose con- 
centrations with respect to time) are statistically signifi- 
cant (here, with respect to different temperatures of the 
medium and preconditioning) or simply due to 



experimental and biological noise or intrinsic variability. 
This objective is generically of great interest in systems 
biology, because it directly affects the construction of 
mechanistic models underlying the investigated time 
courses [4]. In particular, the statistical outcome of the 
analysis forms the basis for the inclusion or exclusion of 
variables in such models, and statistical significance will 
directly lead to new hypotheses as to what biological 
factors or mechanisms might be causing the differences. 

In the illustrative example here, the secondary specific 
task is to determine whether environmental (heat) pre- 
conditioning during population growth affects the 
dynamics of glucose uptake later in the life of the yeast 
cultures. As expected, glucose consumption increases 
with the amount of biomass. However, an unanswered 
question is whether the uptake profile is significantly 
different under altered environmental conditions. Thus, 
we intend to analyze the following questions: First, if 
one adjusts for the differences in initial biomass, are the 
glucose uptake profiles essentially the same or do they 
significantly differ under optimal, heat stress and recov- 
ery conditions? And secondly, do the uptake profiles dif- 
fer between cells that were either grown under optimal 
conditions (30°C) or preconditioned by heat stress (39° 
C) during growth? These questions are not trivial to 
answer, because standard methods of statistics are not 
usually designed to test hypotheses regarding entire 
time profiles for data available at only few time points. 
Formulation of specific hypotheses 

The first specific hypothesis of interest is that the 
dynamic profile of glucose consumption depends on the 
current ambient temperature. The second hypothesis is 
that glucose uptake profiles for cells grown under opti- 
mal conditions differ from the cells grown under heat 
stress conditions (preconditioning). In order to test 
these hypotheses, we need to select a stochastic model 
that fits the dynamics of the rate of change of glucose 
concentration with respect to time. Effects caused by 
moderate differences in initial biomass are accounted 
for by formulating this factor as a covariate in the sto- 
chastic models. 

To facilitate the selection of a suitable stochastic 
model for the glucose uptake dynamics, we first plot the 
rate of change of glucose concentration (mM/min) with 
respect to centered time, i.e., the difference between the 
actual time recorded and the mean time of the experi- 
ment (in minutes), under the three different tempera- 
tures 30°C, 39°C, 30°C for all experiments (Figure 2). 
These plots indeed seem to suggest that initial biomass 
(see dataset Dl) as well as both the ambient tempera- 
ture and the preconditioning with heat during growth 
might have an effect on the glucose uptake dynamics. 
Specifically, it appears that heat stress is associated with 
a higher rate of change in glucose concentration. 
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Figure 1 Glucose uptake profiles for the six experiments under optimal temperature, heat stress, and recovery temperature 



Although there are considerable variations in the 
uptake curves among experiment types, essentially all 
the plots follow a reverse S-shape if the first one or two 
data points are ignored. These first points correspond to 
a delay of uptake initiation that is caused by the experi- 
mental set-up and the fact that the cell population has 
to adjust to the new conditions. In the following, the 
early data points will be included in the analysis, but the 
initial rise in the rate of glucose uptake will not be con- 
sidered as structurally important. 
Model selection 

It is evident from the plots of the data that a nonlinear 
model is required to capture the dynamics of the time 
courses. Accommodating this requirement, as well as 
statistical feasibility, we chose nonlinear mixed-effects 
(NLME) models for our analysis, which are best known 



for analyses of repeated measures and, in particular, 
longitudinal data. In NLME analysis, some or all of the 
fixed and random effects occur nonlinearly in the model 
function. NLME models can be regarded as extensions 
of linear mixed-effects models in which the conditional 
expectation of the response, given the fixed and random 
effects, is allowed to be a nonlinear function of the coef- 
ficients. At the same time, NLME models can also be 
viewed as extensions of nonlinear regression models for 
independent data, in which random effects are incorpo- 
rated in the coefficients to allow them to vary by experi- 
ment type (hereafter referred to as type), thus inducing 
correlation within the types. The latter aspect is particu- 
larly pertinent for our modeling task. 

The NLME model for repeated measures proposed by 
[8] can be thought of as a hierarchical model, where the 
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Figure 2 Rate of change of glucose concentration (mM/min) with respect to centered time (in minutes); under the three different 
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/ observation in the i th type and corresponding to fc* 
level of a covariate is modeled as 



Yijk = f[<Pij, Vijh) + stjk, i = 1, 2, M, 



) ■ 



1,2, 



, rii, k = 1, 2, /. 



(1.1) 



Here M is the number of experimental types, « ; is the 
number of observations on the i th type, / describes the 
levels of the primary covariate, / is a general, real-valued, 
differentiable function with a type-specific parameter 
vector 4>ij and a covariate vector Vj,*, and is a nor- 
mally distributed within-type error term. The function / 
is nonlinear in at least one component of the type-speci- 
fic parameter vector fai, which is modeled as 



cp^Aijp+Bijh b,-~N(0,VO. 



(1.2) 



In this formulation, P is a ^-dimensional vector of 
fixed effects and b t is a ^-dimensional random effects 
vector associated with the i type with variance-covar- 
iance matrix The matrices Ay and By are of appropri- 
ate dimensions and depend on the type and possibly on 
the values of some covariates at the /' observation. 

The model (1.2) was first proposed in [9]. It assumes 
that observations corresponding to different experimen- 
tal types are independent and that the within-type errors 
Zyjc are independently distributed as N(0,a 2 ) and inde- 
pendent of the bi. However, if necessary, this assump- 
tion of independence and homoscedasticity for the 
within-type errors can be relaxed. 

In accordance with the general shape of the time 
course data, we postulate three possible nonlinear mod- 
els for f, namely a simple exponential decay curve, a 
three-parameter logistic model, and a four-parameter 
logistic model. 
Notation 

As in Eq. (1.1), let i denote the index for the experiment 
type, and j the index for the observation number within 
each experiment type. Thus i = 1,2,3,4,5,6 represents the 
six experiment types and /' = 1,2,..., w,-, where, n t is the 
number of observations in the i th experiment type. 
Furthermore, let j,y denote the rate of change of glucose 
concentration for experiment type i of the j th measure- 
ment and tij the centered time recorded for the i* 
experiment type and the j th measurement. In the mathe- 
matical representation below, as well as in the model 
building and selection process, we initially ignore the 
effects due to temperature and preconditioning during 
cell growth. However, the effect of these covariates, 
namely temperature, preconditioning and the tempera- 
ture-preconditioning interaction, is incorporated in the 
best model and presented later (see Eq. 3.5). 
Model 1: Exponential Decay Curve The exponential 
model has plenty of theoretical and empirical support. 



In particular, if the dynamics of glucose (y) followed a 
simple transport process (from the medium into the 
cells) or a first-order kinetic process, it would be 
described as a linear differential equation of the type 
y = —k.y, whose solution and rate of change are expo- 
nential functions. Although, the exponential decay 
model is common for the first-order kinetic equations, 
it can fail outside the range of the data. Hence we chose 
instead a more flexible model with more theoretical and 
empirical support, namely a three-parameter exponential 
decay curve that does not necessarily have zero as its 
asymptote. Thus, in this first case we express the rate of 
change of glucose concentration yu of experiment type i 
at time t t j for i = 1,2,3,4,5,6 / = 1,2,..., «; 



YH 



(2.1) 



The function contains three parameters, of which, Xi 
is the right asymptote, X 2 is the total amount of glucose 
ultimately taken up, and X 3 is the time needed to take 
up half of the maximum amount [10]. 8,y is a normally 
distributed {N{0,a )) error term, which is assumed to be 
identically and independently distributed (i.i.d.). Xi and 
X 2 are linear parameters while X 3 is a nonlinear para- 
meter. Additional file 1 Table SI in the Appendix, 
shows that the three-parameter model fits the data sig- 
nificantly better than the simple two-parameter model. 
Model 2: Three-parameter Logistic Model The three- 
parameter logistic model expresses the rate of change in 
glucose concentration yu of experiment type i at time t,y as 



YH = M/(l 



-ftr**)As) + £ ... 



(3.1) 



Like the previous one, this model is governed by three 
parameters. X 1 is the maximum rate of change for glu- 
cose concentration X 2 is the time at which the maxi- 
mum rate of change for glucose concentration attains 
half of its maximum value, and X 3 is the reciprocal of 
the decay rate. (As generic support for the appropriate- 
ness of this model, see [9] and [10]). As before, e 2 y is an 
i.i.d. N(0,a 2 ) error term. X l is a linear parameter, while 
X 2 and X 3 are nonlinear parameters. 

The format in Eq. (3.1) is not often seen in biochemis- 
try, where nonlinear models are preferred if they can be 
interpreted mechanistically and if their parameters have a 
natural physical interpretation. For instance, a typical 
description of the glucose uptake process might be in the 
form of a kinetic Michaelis-Menten process of the type 



dy Vmaxy 

dt K M +y 



(3.2a) 



This differential equation can be solved analytically, 
but the explicit solution of y as a function of time is 
rather convoluted. However, the simpler three- 
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parameter logistic model (3.1) offers a very good 
approximation. For instance, if the parameters of the 
Michaelis-Menten function (3.2*) are arbitrarily chosen 
as V max = 2.4 K M = 3.3 and y{0) = 8, the logistic model 
(3.1) with X 1 = 11.6, X 2 = 1-4 and X 3 = -1.85 results in a 
graph that is essentially indistinguishable from the solu- 
tion of the differential equation in (3.2*) (Figure 3). 
Thus, within minimal experimental noise, the simple 
three-parameter logistic function models a Michaelis- 
Menten-type uptake function with sufficient accuracy. 
Model 3: Four-parameter Logistic Model The four- 
parameter logistic model is given by the mathematical 
formula 

Yij = X 1 + [X 2 - X 1 )/{1 + e CXs-%)/*.) + e ., (4.1) 

Here, is the horizontal right asymptote, X 2 is the 
horizontal left asymptote (baseline) and X 3 is the time 
that corresponds to the rate of glucose consumption 
halfway between the asymptotes ([9] and [10]). X 4 is the 
reciprocal of the decay rate and £y is again an i.i.d. N(0, 
a ) error term. In Eq. (4.1), Xi and X 2 are linear para- 
meters, while X 3 and X 4 are nonlinear parameters. To 
resolve the issue of lack of identifiability, we require that 
X 4 must be strictly positive. Clearly, (3.1) is a special 
case of (4.1) with X 1 = 0. 

The parameters of the models (2.1), (3.1), and (4.1) are 
obtained with nonlinear least squares (nls) regression. 
To fit a separate nonlinear model to each dataset and 
thus allowing the individual effects (experiment types 
including temperatures) to be incorporated in the para- 
meter estimates, we express the models (2.1), (3.1) and 
(4.1) as 



yy = X n + kae-**'** + sy. (2.2) 
ro-Afl/fl + e-Gr^Aaj + efl. (3.2) 

y tj = X n + [X i2 - ka)/(l + + ey. (4.2) 

where, as before, the e,y's are i.i.d. N(0,a 2 ) errors. The 
functions were subsequently used in our NLME model 
with a mixture of fixed and random effects. In the statis- 
tical software package R [11], models (2.2), (3.2), and 
(4.2) are considered nlsList models, since their para- 
meters are estimated with nonlinear least squares 
methods. 

To account for fixed and random effects in a non- 
linear mixed-effects model, it is useful to reformulate 
the models (2.2), (3.2) and (4.2) respectively as 

y {j = Xi + (An - 

[X 2 + [X i2 - l 2 )\ e -^+^-^)\ + ey' { ' 

Yij = &i + (X a -Xi))/ 

(1 + e-(%-[i2+(Ae-*a)]/[i3+(AB-*3)l)) + Sij' 

ytj = (Xi + {Xn - Xi)) + ([X 2 + 

[X a - X 2 )] - [X! + [Xa - , (4.3) 

where the x quantities denote the means of the para- 
meters for each type. In nlsList model(s), the deviations 
of the coefficients from their means are treated as para- 
meters that need to be estimated. Mixed-effects models, 
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on the other hand, represent these deviations from the 
mean value as random effects, treating each of the types 
as a sample from a population. The nonlinear mixed- 
effects versions of (2.2), (3.2) and (4.2) are therefore 

Yij = (Pi + ba) + (Pi + b t2 )e-"^ ba) + s l} 



(2.4) 
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Yij = (Pi +hi) + ((Pi + b i2 ) - (Pi + ba.))/ 

(1 + e ((fe+fcs)-%)/(A+*i4))^ + g.. 



vw 



'/ 

12 
i3 



(4.4) 



The fixed effects Pi,P2>P3 (and PJ represent the mean 
values of the parameters for each type, X, is the popula- 
tion of experiment types and the random effects, bi\,ba, 
b i3 (and b^), represent the deviations of the X { from 
their population average. X/s contain both fixed and 
random effects. The random effects are assumed to be 
distributed normally with mean 0 and variance-covar- 
iance matrix y. Random effects corresponding to differ- 
ent types are assumed to be independent. The within- 
type errors s,, are assumed to be independently distribu- 
ted as N(0,a ) for different i, j and to be independent of 
the random effects. Models (2.4), (3.4), and (4.4) are 
considered nonlinear mixed-effects (NLME) models 
since both fixed and random effects are incorporated in 
the models. The nonlinear mixed-effects models (2.4), 
(3.4), and (4.4) offer a compromise between the rigid 
nonlinear least squares (nls) models (2.1), (3.1) and (4.1) 
and the over-parameterized nlsList models (2.2), (3.2) 
and (4.2). They accommodate variations for each experi- 
mental type through the random effects, but tie the dif- 
ferent types together through the fixed effects and the 
variance-covariance matrix \|/. A crucial step in the 
model-building of mixed-effects models is deciding 
which of the coefficients (parameters) in the model need 
random effects to account for their between-type varia- 
tion and which can be treated as purely fixed effects. 



We now fit a nonlinear regression model for the rate 
of change of glucose concentration with respect to time 
for each temperature. For this purpose, we require for 
each temperature the same starting time for glucose 
uptake. Since we are interested in modeling and drawing 
statistical inferences for the rate of change in glucose 
concentration and the collection times are equally 
spaced for each temperature, we consider the starting 
times for each temperature to be the ones for which the 
first glucose concentration data are available. Further- 
more, it simplifies the computation, without loss of gen- 
erality, if we use centered times, by considering the 
differences between recorded times and the mean time 
of the experiment type as in Figure 2. 

Initially ignoring the effects of temperature, we fit the 
nonlinear least squares model of the data to determine 
which of the three model functions ((2.1), (3.1) or (4.1)) 
best fits the time trends in the data. Then we fit non- 
linear least squares curves for each experiment type in 
order to determine if the models require random effects 
(nlsList models (2.2), (3.2), (4.2)). The main purpose of 
this preliminary analysis provided by nlsList is to suggest 
a structure for the random effects to be used in a non- 
linear mixed-effects regression model. Another impor- 
tant advantage of using an nlsList object is that it 
automatically provides initial estimates for the fixed 
effects, the random effects, and the random-effects cov- 
ariance matrix. 

Following the model-building process for nonlinear 
mixed-effects models proposed by [9], we first checked 
all assumptions regarding the random effects, tested for 
heteroscedasticity of the within-type error variance and 
the correlation of the random effects. We then per- 
formed a likelihood ratio test, and checked the Akaike 
Information Criterion (AIC) as well as the Bayesian 
Information Criterion (BIC) to select the nonlinear 
regression model that best fits the data. In this case, it 
turned out to be the three-parameter logistic model. 
The model parameters were obtained with the method 
of maximum likelihood estimation, and the parameter 
estimates are reported by fitting the best nonlinear 
regression model with fixed effects terms; temperature, 
preconditioning and the temperature-preconditioning 
interaction included in the model as covariates. Thus, 
the mathematical formula for the three-parameter logis- 
tic nonlinear model with temperature, preconditioning 
and temperature-preconditioning interaction term as 
covariates for all the parameters is given by 



Yij = (Pi + Yn x m + YuXm + e x x i2 + 

%XlXi\lXi2 + %12XiiiXi 2 )/(l + 

-C%-CA+V4l»iU+K!2*<12+S2«i2+f21«ill'Ci2+fe!2*jl25Cj2))/ 
0O33+ra*ll+W2*<12+03*i2+fel*m*a+l32*fl2*i2) J + Ej - 



(3.5) 
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11 Tempi = heat 
0 Tempi = re cov ery 

{1 Tempo = re cov ery . 
rj i and 
0 Temp 2 = heat 

{1 precond = present 
0 precond = absent 

where, as before,Pi, p 2 and p 3 represent the mean 
values of the parameters for each type Xi in the popula- 
tion of experiment types and the random effects, and b a , 
ba and ba represent the deviations of the X, from their 
population average. Furthermore, let k denote the index 
for the covariate temperature used while observing glu- 
cose consumption. Thus, Temp k represents the specific 
temperature index k = 1,2 that corresponds to heat stress 
or recovery conditions, respectively. Note that the tem- 
perature 'optimal' is used as reference category. The 
index precond represents the two preconditioning levels 
corresponding to either optimal or heat stress conditions 
during the late growth of the cell populations, before the 
actual experiments started, yn, y 2 i and 731 (712, Y22 and 
y 32 ) represent the coefficients for the difference between 
the optimal temperature category and heat stress (or 
recovery). 81, 82 and 83 represent the coefficients for the 
difference between optimal condition and precondition- 
ing during initial growth. £ U) § 2 i and (^ 12 , ^22 and ^ 32 ) 
represent the coefficients for the interaction between the 
heat stress-preconditioning and the temperatures while 
observing glucose consumption. 

Methods for Testing the Temperature Effect 

Analysis of variance (ANOVA) was used to test if tem- 
perature has an effect on the glucose dynamics. For 
detecting the overall difference in the dynamics due to 
differences in temperatures, we considered two models, 
namely a full model, in which parameters of the non- 
linear regression model vary with respect to temperature 
and biomass concentration, and a reduced model with 
constant and temperature independent parameters. The 
likelihood ratio test was performed for assessing the 
overall temperature effect, and characterized with a p- 
value. 

For pair-wise comparisons, we considered appropriate 
subsets of the data. Initially, we used ANOVA with the 
approximate F-test to test the joint significance of the 
fixed effects introduced in the model. The F-test for 
joint significance was performed, and j?-values were 
reported for the test as well. 

When using the F-test, the assumption of normality 
should be justified. However, it is extremely difficult to 



check for normality in the case of nonlinear mixed- 
effects models. In fact, as the random effects are them- 
selves not observable, a direct check of their normality 
is not possible. Moreover, the complicated way in which 
the random effects enter a nonlinear mixed model 
makes it impossible to check the assumption of normal- 
ity of the errors. An interesting discussion on this speci- 
fic topic can be found in Hartford and Davidian (see 
[12]). These authors mention the fact that, because the 
random effects are unobservable latent model compo- 
nents, no straightforward diagnostic is available to check 
the validity of the normality assumption. However, huge 
simulations have shown that estimation of fixed para- 
meters may not be severely compromised (e.g., [13]) 
even if the assumption of normality of the random 
effects is violated. 

Methods for Testing the Effect of Preconditioning 

To investigate the effect of heat preconditioning during 
growth, we used a similar procedure as before for test- 
ing the effect of temperature on glucose uptake 
dynamics, namely the ANOVA likelihood ratio test. In 
addition, the approximate F-test was used for testing the 
joint significance of fixed effects terms in the model. 

Simulation Study 

To illustrate and justify the overall usefulness of the 
proposed method of fitting a three-parameter logistic 
model to the glucose uptake profiles we carried out a 
robustness analysis with the help of a simulation study. 
The steps in the simulation study were as follows: 

In the first step, we simulated £*■ ~ N(0, a 2 ) i = 
1,2,3,4,5,6; / = 1,2,. We used the estimated value of 
a = 0.3103 from the dataset. We further simulated the 
random effects b* v b* 2 , b* 3 ~ N(0, i/). We used multiples 
of the estimates of the variance-covariance matrix of the 
random effects as obtained from the real data as \\i to 
generate the simulated data. We also used the parameter 
estimates of the fixed effects coefficients directly from 
the data to generate the simulated data. Thus, the simu- 
lated data [y*j, Uj) were generated from the model 

y * = (/Jj + + el-to-Cft^JVCft^)]) + e ... 

A sample size of 1000 was used. The parameters were 
estimated from each simulated data set. We systemati- 
cally changed the random effects variances simulated 
data sets multiple times, and estimated the parameters 
of the three-parameter logistic model for every instance. 
Specifically, parameter estimates of the fixed effects 
terms, bias and standard error for different choices of 
\|/s are reported. 
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Results 

Results for Model Fitting and Statistical Methods of 
Significance Analysis 
Model Fitting 

We employed all NLME models with the three-para- 
meter exponential and the two logistic functions to fit 
all data simultaneously. Among them, the three-para- 
meter logistic model had the lowest AIC and BIC values 
of 1911.637 and 1953.981, respectively (see Additional 
file 1 Table S2, Additional file 1 Table S3, and Addi- 
tional file 1 Table S4 in the Appendix). All but one of 
the parameter estimates (intercept not included) of this 
model were found to be significant at the 5% signifi- 
cance level (Table 1), indicating that the three-para- 
meter logistic model provides overall the best 
simultaneous fit for all data. Indeed, the logistic mixed- 
effects nonlinear regression model fits the glucose 
uptake dynamics through the rate of change of glucose 
concentration rather well for all given datasets (Figure 
4). The (initial) temperature condition 'optimal' is con- 
sidered as the reference temperature and referred to as 
the intercept of the model. 

The intercept terms (Pi. Intercept, p 2 . Intercept, and p 3 . 
Intercept) represent the mean parameter estimates for 
the optimal temperature condition for the six experi- 
ment types. Specifically, pi. Intercept is the maximum 
value for the rate of change of glucose concentration for 
the optimal temperature averaged over the six experi- 
ment types. For the optimal temperature, P2.1ntercept is 
the time at which the uptake rate curve attains half of 
its maximum value averaged over the six experiment 
types. The estimated P 3. Intercept parameter represents 
the average rate of glucose utilization in the six experi- 
ment types for the optimal temperature. 

The parameter estimates associated with an increase 
in temperature from optimal to heat stress for the six 
experiment types are quantified by P^Heat, p 2 .Heat and 



Table 1 Parameter estimates for the three-parameter 
logistic model with temperature included in the model 
as covariate. 



Parameter 


Estimate 


t-value 


p-value 


P,. (Intercept) 


1 0.4096 


16.1326 


0.0000 


P,.Heat 


0.6456 


1.1436 


0.2533 


Pi. Recovery 


-2.0912 


-4.1018 


0.0000 


p 2 .(lntercept) 


-0.8622 


-1.5464 


0.1227 


p 2 .Heat 


-0.9627 


-6.9235 


0.0000 


p 2 .Recovery 


-0.7040 


-3.2702 


0.001 1 


p 3 . (Intercept) 


-1 .0846 


-6.0598 


0.0000 


p 3 .Heat 


0.5000 


6.4934 


0.0000 


P 3 .Recovery 


-0.3239 


-2.6566 


0.0081 



The optimal condition is considered as the reference temperature. 



P 3 .Heat. The parameter estimate Pi.Heat, on average for 
the six experiment types, accounts for an increase of 
0.65 in the maximum rate of change in glucose concen- 
tration under heat stress compared to optimal tempera- 
ture, which is not significant at the 5% significance level 
{p = 0.2533). However, the p 2 .Heat parameter estimate 
is negative and larger than p 2 . Intercept, thereby imply- 
ing that on average the time at which the rate curve 
attains half of its maximum due to heat stress is shorter 
than under optimal temperature conditions. The esti- 
mate is negative with a very low p-value (p < 10" 4 ), 
which implies that glucose is taken up significantly fas- 
ter during heat stress compared to the optimal tempera- 
ture. The time interval remains the same for the 
optimal and the heat stress category, which indicates 
that the decay rate of the three parameter logistic curve 
is larger in the heat stress category. This conclusion is 
supported by the estimate of parameter p 3 .Heat, which 
is positive and significant in heat stress with a very 
small p-vahie. This result implies that the rate of decay 
during heat stress is faster than under optimal tempera- 
ture conditions. 

The effect of a return to the optimal temperature dur- 
ing recovery after heat stress is captured by the para- 
meter estimates Pi. Recovery p 2 . Recovery and p 3 . 
Recovery. The estimates are negative in sign with very 
low /7-values implying that the glucose uptake profile 
during recovery has on average a lower maximum value 
as depicted by Pi. Recovery, a shorter time period to 
attain half of its maximum (p 2 . Recovery) and a slower 
rate of decay (p 3 . Recovery). This translates into a flatter 
uptake rate curve during recovery as compared to the 
initial phase of optimal ambient temperature. 

Results of the Significance Analysis 

Table 2 shows the results of an ANOVA for testing the 
effect of ambient temperature on the rate of change of 
glucose consumption over time. The likelihood ratio test 
results in a very small p -value (<0.0001), which again 
suggests with high confidence that temperature signifi- 
cantly affects the glucose uptake dynamics. 

For the second part of the analysis, testing the effect 
of heat preconditioning on the glucose uptake dynamics 
over time, the likelihood ratio test for the ANOVA 
results in a very small p-val\ie (p = 0.0004) indicating 
that preconditioning significantly affects glucose uptake 
dynamics. This result also suggests that the three-para- 
meter logistic mixed-effects model explains the glucose 
uptake characteristics for this situation as well. 

To investigate the overall effect of preconditioning, we 
considered two models as before: a full model and a 
reduced model. In the full model the parameters 
obtained with nonlinear regression depend on precondi- 
tioning, and temperature, whereas in the reduced model 
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Figure 4 Plots of observed data (blue circles; adjusted) and time courses (red lines) obtained from a simultaneously fitted three- 
parameter logistic model in the implementation of a reduced model of Eq. (3.5) with temperature as the only covariate No 

measurements were taken during recovery in experiment NG3. Notice that NG1 (right column), while normally grown, started with a 
substantially higher biomass than NG2-NG4. PC1-PC2 indicate preconditioned cells. 



the parameters depend only on temperature, which is 
accounted for as the primary covariate. The results are 
again significant (Table 3). 

We can furthermore investigate the interaction effect 
between temperature and preconditioning. In this case, 
the full model is formulated in such a fashion that the 
parameters of the model depend on temperature, precon- 
ditioning and the temperature-preconditioning interac- 
tion effect, whereas the parameters of the reduced model 
depend on temperature and preconditioning as the only 
fixed effects terms. Our results indicate that the interac- 
tion effect of temperature and preconditioning is highly 
significant (p < 0.0001) at the 5% significance level, which 
implies that preconditioning and later heat stress interact 
significantly and lead to a combined effect on the glucose 
uptake profiles (see Additional file 1 Table S5 in the 
Appendix for ANOVA results). The three-parameter 
logistic mixed-effects model again explains glucose 
uptake dynamics very well for this situation. 



The effect of initial biomass on glucose uptake 
dynamics was investigated following a similar procedure 
as above. Our results show that the p-velue for the like- 
lihood ratio test is very small (p < 0.0001) indicating 
that the amount of cell biomass at harvest significantly 
affects glucose uptake kinetics (see Additional file 1 
Table S6 in the Appendix for ANOVA results). For this 
scenario, the parameters of the full three-parameter 
logistic mixed-effects nonlinear regression model vary 
with respect to temperature and the amount of biomass 
harvested, while the parameters in the reduced model 
depend only on temperature. 

Pair-wise Comparisons 

In addition to the overall assessments, it is also possible 
to execute pair-wise comparisons between different tem- 
perature conditions during glucose uptake. For these 
comparisons, one condition needs to be selected as refer- 
ence, and as before we choose the optimal temperature 



Table 2 Results of ANOVA and the overall likelihood 
ratio test for determining whether the effect of 
temperature on glucose uptake is significant. 

Model Df AIC BIC logLik Test L.Ratio p- 

value 

Full 36 937.4367 1089.876 -432.7186 



Table 3 Results of ANOVA and the overall likelihood 
ratio test for determining whether heat stress during 
growth significantly affects the glucose uptake dynamics. 

Model Df AIC BIC logLik Test L.Ratio p- 

value 

Full 24 990.3122 1091.938 -471.1561 



Reduced 30 1117.0356 1244.068 -528.5178 1 vs 191.5988 <0.0001 Reduced 21 1002.7044 1091.627 -480.3522 1 vs. 18.3923 0.0004 

2 2 
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for this purpose. The analysis is achieved with an approx- 
imate F-test for testing the joint significance of the terms 
in the models. The main results consist of very small p- 
values (p <0.0001) that indicate significant differences 
between all pairs of conditions (Table 4). The low up- 
values in all three cases demonstrate that the added tem- 
perature terms in the final models are highly significant 
in all cases. Interestingly, the glucose uptake dynamics 
also differs between optimal and recovery conditions. 

The significance of including temperature in the 
model indicates that a substantial part of the type-to- 
type variation in the dynamics of glucose uptake is 
explained by pair-wise differences in temperature. In 
other words, the fixed effects terms introduced in the 
model to explain the temperature variability in the para- 
meter estimates are significantly different from zero at 
the 5% significance level. Since the result of the F-test is 
highly significant, we can employ an additional i-test in 
order to identify exactly which parameters are signifi- 
cant and which are not in explaining the type-to-type 
variability in the model. Table 4 also identifies the p- 
value for the approximate F-test for testing the effect of 
preconditioning in the model as very small (p = 0.0033), 
which indicates that the fixed effect term (optimal or 
heat stress during growth), which was introduced into 
the model to explain type-to-type variability due to pre- 
conditioning, is significantly different from zero. 

As indicated earlier, heat stress conditions during cell 
culture growth affect the dynamics of glucose consump- 
tion. In Table 5 the intercept terms (pi. Intercept, P2-Inter- 
cept, and p 3 .Intercept) represent the parameter estimates 
averaged over the six experiment types under optimal con- 
ditions, whereas the preconditioning terms (p x . Precondi- 
tioning, p 2 . Preconditioning, and p 3 . Preconditioning) 
represent the corresponding parameter estimates for com- 
paring cell cultures grown under heat conditions relative 
to those grown under optimal conditions and averaged 
over the six experiment types. For instance, the parameter 
estimate p x . Preconditioning indicates the difference in the 
maximum value for the glucose uptake rate curve averaged 



Table 4 Approximate F-test; for joint significance of the 
fixed effects terms (temperature and preconditioning) in 
the model. 



Fixed Effect 


Comparison 


NumDF 


denDF 


F-value 


P- 

value 


Temperature 


Optimal & Heat 
Stress 


3 


325 


1 1 8.6734 


<0.0001 


Temperature 


Optimal & 
Recovery 


3 


331 


51.3787 


<0.0001 


Temperature 


Heat Stress & 
Recovery 


3 


331 


56.1397 


<0.0001 


Preconditioning 


Absence & 
Presence 


3 


499 


4.6458 


0.0033 



Table 5 Parameter estimates for the three-parameter 
logistic model with preconditioning included in the 
model as covariate. 



Parameter 


Estimate 


t-value 


p-value 


p -^Intercept) 


1 1 .5740 


11.3163 


0.0000 


p q .Preconditioning 


-3.0639 


-1.7608 


0.0789 


p 2 .(lntercept) 


-2.2904 


-3.8724 


0.0001 


p 2 .Preconditioning 


24782 


2.3383 


0.0198 


p 3 .(lntercept) 


-0.8490 


-5.1658 


0.0000 


P 3 .Preconditioning 


-0.6293 


-1.9784 


0.0484 



Cell cultures grown under cold conditions are considered as the reference 
category. 



over the six experimental types for cultures grown under 
heat conditions compared to cells grown under optimal 
temperature conditions. Its negative value indicates that 
the average maximum for the rate of glucose utilization is 
lower for preconditioned cells than for control cells grown 
under optimal condition, implying a lower uptake rate for 
cultures grown under heat conditions compared to those 
grown under optimal conditions. The parameter estimate 
P 2 . Preconditioning reflects the difference in the time at 
which the maximum rate of glucose utilization attains half 
of its value. Again, this value results from averaging over 
the experimental types for the preconditioned cells com- 
pared to controls. Its positive value indicates that, on aver- 
age, the preconditioned cells take more time to consume 
glucose and hence had a slower glucose uptake than the 
control cells, which confirms a result from an earlier mod- 
eling study that used entirely different methods of analysis 
[4]. The parameter estimate p 3 . Preconditioning is the 
mean difference in the scale of glucose uptake for the cells 
grown under heat versus control conditions. 

Results of the Simulation Study 

In order to study the robustness of the parameter esti- 
mates using the nonlinear mixed-effects model for the 
rate of glucose uptake that was deemed best (the three- 
parameter logistic model) we used a simulation study. 
We checked the performance of the results over multi- 
ple variance-covariance matrices y and report the 
results for two representative cases where the variance- 
covariance parameter matrices are multiples of y, which 
were estimated directly from the data: 

/0.001159 0 0 \ 

ir 2 = 0.001 * f = 0 0.001313 0 

\ 0 0 0.000398/ 

and 

/0.011588 0 0 \ 

f 3 = 0.01 * if = 0 0.013126 0 

\ 0 0 0.003982/ 
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Table 6 Parameter estimates for the simulation study 
with the three parameter logistic model with random 
effect generated from N (0, 0.001 "estimated \\i from the 
data). 



Parameter 


Simulation 
parameter 
values 


Simulated 
parameter 
estimate 


Estimated 
Bias 


Estimated 

Standard 

Error 


Pi.lntercept 


10.4096 


10.4321 


-0.0225 


0.0724 


Pi. Heat 


0.6456 


0.6146 


0.0311 


0.0956 


Pi. Recovery 


-2.0912 


-2.1051 


0.0138 


0.1058 


p 2 .lntercept 


-0.8622 


-0.8665 


0.0043 


0.0347 


p 2 .Heat 


-0.9627 


-0.9559 


-0.0068 


0.0399 


p 2 .Recovery 


-0.7040 


-0.6986 


-0.0055 


0.0663 


p 3 .lntercept 


-1.0846 


-1.1007 


0.0161 


0.0277 


p 3 Heat 


0.5000 


0.5128 


-0.0128 


0.0320 


p 3 . Recovery 


-0.3239 


-0.3332 


0.0093 


0.0524 



The estimated biases and standard deviations of the 
fixed effects parameters are reported in Tables 6 and 7 
respectively. Estimated biases and standard deviations 
both remain quite small and are therefore acceptable. 
For other values of the parameters, the parameter esti- 
mates are also robust, which suggests that the method 
works reasonably well (results not reported). 

Discussion and Conclusions 

We have developed stochastic non-linear regression 
models to fit and analyze sparse biological time course 
data. In the application shown here, these data were 
generated with NMR methods. They describe the 
dynamics of glucose uptake by yeast cells at different 
temperatures. Some of the cells in the study were grown 
under optimal control conditions, while others were pre- 
conditioned with heat treatment during growth. To 
facilitate strict comparisons, we used centered time, 

Table 7 Parameter estimates for the simulation study 
with the three parameter logistic model with random 
effect generated from N (0, 0.01 *estimated \\i from the 
data). 



Parameter 


Simulation 
parameter 
values 


Simulated 
estimate 


Estimated 
Bias 


Estimated 

Standard 

Error 


Pi.lntercept 


10.4096 


10.3636 


0.0460 


0.0595 


B v Heat 


0.6456 


0.6663 


-0.0207 


0.0802 


P,. Recovery 


-2.0912 


-2.0240 


-0.0673 


0.0853 


p 2 .lntercept 


-0.8622 


-0.8440 


-0.0182 


0.0275 


P 2 .Heat 


-0.9627 


-0.9732 


0.0105 


0.0320 


p 2 .Recovery 


-0.7040 


-0.6898 


-0.0143 


0.0512 


P 3 .lntercept 


-1.0846 


-1.1172 


0.0326 


0.0224 


P 3 Heat 


0.5000 


0.5040 


-0.0039 


0.0246 


P 3 .Recovery 


-0.3239 


-0.3288 


0.0050 


0.0387 



which took advantage of the fact that the time intervals 
for data collection were the same for all experimental 
set-ups and avoided problems caused by different start- 
ing times in the various experiments. 

The main goal of the study was to develop rigorous 
statistical tests to assess the effects of temperature and 
preconditioning on the observed temporal glucose 
uptake profiles. To achieve this goal, we designed non- 
linear mixed-effects regression models and analyzed 
them with customized ANOVA and maximum likeli- 
hood ratio tests. 

The results indicate that a change in temperature 
from optimal to heat stress conditions (30°C to 39°C) 
produced significant differences in glucose uptake pro- 
files. Specifically, the increase in temperature had a posi- 
tive effect [14] on the rate of glucose uptake, while a 
corresponding decrease in temperature from heat stress 
to recovery resulted in a reduction in rate of glucose 
uptake. However, the recovery profiles were found to be 
different from the initial profiles (initial optimal condi- 
tion), although the temperature was the same. We also 
determined that preconditioning with heat during cell 
growth resulted in significant differences in the glucose 
uptake profiles later in life. These differences presum- 
ably reflect a preparation of the cells to survive similar 
stresses [15,16] and [17] and to recover in a timely man- 
ner. From a methodological point of view, the results 
indicate that the methods described here are able to 
detect subtle difference in time course data from normal 
and stressed cells. 

The statistical methods presented in this paper reach 
far beyond tests for the effects of different stress condi- 
tions. In fact, it appears that the exact same statistical 
procedures should be applicable and very useful for 
finding significant differences in any short biological 
time course data of a similar format. For instance, it 
seems that temporal expression profiles of genes or pro- 
teins from high-throughput experiments could be clus- 
tered by pathways and then analyzed with the proposed 
methods to discover differences between normal and 
disease conditions. 

The estimation of parameter values from biological 
time series data is not a trivial problem, and the estima- 
tion algorithms for nonlinear mixed-effects are compu- 
tationally complex and often provide less accurate 
inference than for linear effects models. Furthermore, 
nonlinear regression models require starting estimates 
for the fixed effects coefficients. Determining reasonable 
starting parameter estimates is somewhat of an art and 
not intuitive in many situations. Nevertheless, the con- 
cepts presented here are quite general, and algorithms 
exist in statistical software packages like R for their 
implementation and analysis. Finally, this type of analy- 
sis, like almost all dynamic modeling efforts in systems 
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biology, requires the choice of a parametric model. It is 
unclear how one could get around this requirement, 
which is a challenge for about any modelling study in 
systems biology. One could potentially use canonical 
forms, such as piecewise power-law functions [18], but 
while these would be slightly more generic, they would 
still be parametric. 

Additional material 



Additional file 1: Appendix. This file provides additional model 
comparisons tables. These tables explain the choice of best model, the 
test of interaction effect, the test of biomass concentration as well as 
comparing the two and three parameter exponential models. 
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